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Abstract 

Closed form option pricing formulae explaining skew and smile are 
obtained within a parsimonious non-Gaussian framework. We extend 
the non-Gaussian option pricing model of L. Borland (Quantitative 
Finance, 2, 415-431, 2002) to include volatility-stock correlations con- 
sistent with the leverage effect. A generalized Black-Scholes partial 
differential equation for this model is obtained, together with closed- 
form approximate solutions for the fair price of a European call op- 
tion. In certain limits, the standard Black-Scholes model is recovered, 
as is the Constant Elasticity of Variance (CEV) model of Cox and 
Ross. Alternative methods of solution to that model are thereby also 
discussed. The model parameters are partially fit from empirical ob- 
servations of the distribution of the underlying. The option pricing 
model then predicts European call prices which fit well to empirical 
market data over several maturities. 
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1 Introduction 



Over the past three decades, since the seminal works of Black, Scholes and 
Merton [1, 2], the basic ideas of arbitrage-free option pricing have become 
quite well-known. The Black and Scholes model is such a standard bench- 
mark for calculating option prices, that is used by traders world-wide in 
spite of the fact that the prices it yields do no match ones observed on the 
market place. This difference is clearly a signature of the fact that, after 
all, the Black-Scholes formula is just a model of reality, based on a set of 
assumptions which obviously only approximate (often quite remotely) the 
true dynamics of financial markets. Most importantly, real markets are in- 
complete and risk cannot in general be fully hedged away, as it is the case 
in the Black-Scholes world [3]. Nevertheless, prices are quoted in terms of 
the Black-Scholes model, and in order to obtain the market prices from this 
model, it is commonplace to adjust the volatility parameter a which enters 
the model. In other words, traders use a different value of a for each value 
of the option strike price as well as for each value of the option expira- 
tion time T. This is tantamount to having a continuum of different models 
(corresponding to each value of a) for the different options on the very same 
underlying instrument. A plot of a over strike and maturity then yields a 
surface, often convex and sloping, which is referred to as the smile or the 
skew, or more generally the smile or skew surface. Given a smile surface, one 
can plug those values of a into a Black-Scholes model and obtain prices for 
each strike and time to expiration. 

Many attempts have been made to modify or extend the Black-Scholes 
model in order to accommodate for the skew observed on option markets. 
A large class of those models is based on modeling the skew surface itself. 
For example, so-called local volatihty models [4, 5] aim to fit the observed 
skew surface by calibrating a function (Tioc{S,t) (where S is the stock price 
and t the time) such that it reproduces actual market prices for each K 
and T. This function aioc{S, t) is then used in conjunction with the Black- 
Scholes model to perform other important operations such as the pricing of 
exotic options, hedging, and so on. The problem is however, that while aioc 
by construction allows one to reproduce the set of prices to which it was 
calibrated, it docs not contain any information about the true dynamics of 
the underlying asset; therefore it can actually lead to worse hedging strategics 
and erroneous exotic prices than would have been obtained with simply using 
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Black-Scholes without any attempts to account for the smile [6]. Nonetheless, 
local volatility models have been widely used by many banks and trading 
desks. 

Levy processes [7, 8, 9, 10], stochastic volatility models [11, 12, 13, 6, 14, 
15] or cumulant expansions around the Black-Scholes case [16, 17, 18, 19, 10] 
constitute approaches which have been successful in capturing some of the 
features of real option prices. For example the recent 'stochastic alpha, beta, 
rho' (SABR) model ([6], and see Appendix C) can be well- fit to empirical 
skew surfaces. It also provides a better model of the dynamics of the smile 
over time. Nevertheless, in all these models the focus is still on fitting or 
somehow calibrating the parameters of the model to match observed option 
prices. The difference in our current approach will instead be to introduce a 
stock price model capturing some important features of the empirical distri- 
bution of stock returns. Our model will be intrinsically more parsimonious 
than stochastic volatility models in the sense that we have just one source of 
randomness, which also allows us to remain within the framework of complete 
markets. The option pricing methodology yields closed form approximate for- 
mulae based on such a model, thereby predicting option prices rather than 
fitting parameters to match observed market prices. We find good agree- 
ment between theoretical and traded prices, lending support to the possible 
validity and potential applicability of the model. 



2 Non-Gaussian Stock Price Model with Skew 

The standard Black-Scholes stock price model reads 



where doo represents a zero mean Brownian random noise correlated in time 
t as 



Here, n represents the rate of return and a the volatility of log stock returns. 
This model implies that stock returns follow a log-normal distribution, which 
is only a very rough approximate description of the actual situation. In fact, 
real returns have strong power-law tails which are not at all accounted for 
within the standard theory. Furthermore, there is a skew in the distribution 



dS — iiSdt -\- aSduj 



(1) 



{duj{t)duj{t')) F = dtdt'5{t - t'), 



(2) 
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such that there is a higher probabihty of large negative returns than large 
positive ones. Both the tails and the skew depend on the time over which 
the returns are calculated. In general, it is observed that the power-law 
statistics of the distributions are very stable, exhibiting tails decaying as —3 
in the cumulative distribution for returns taken over time-scales ranging from 
minutes to weeks, only slowly converging to Gaussian statistics for very long 
time-scales [20, 10]. Similarly, the skew of the distribution varies with the 
time-scale of returns such that it is largest for intermediate time-scales [10]. 

The deviations of the statistics of real returns to those of the log-normal 
model of Eq. (1) become particularly important when it comes to calculating 
the fair price of options on the underlying stock. For example, the simplest 
'European' call option, which is the right to buy the stock S at the strike price 
X at a given expiration time T. The call will expire worthless if S{T) < K, 
and profitably otherwise. The fair price of the call thus depends on the 
probability that the stock price S(T) exceeds and if one uses the wrong 
statistics in the model then the theoretical price will differ quite a bit from 
the empirically traded price. (Interestingly enough, market players intuitively 
adjust the price of options to be much more consistent with the true statistics 
of stock returns [19], even though the log-normal Black-Scholes model is 
widely used by traders themselves to get an estimate of what the price should 
be.) 

In this paper, we develop a model of the underlying stock which is con- 
sistent with both fat-tails and skewness in the returns distribution. We then 
derive option prices for this model, obtaining (approximate) closed-form so- 
lutions for European call options. 

Our model bases on the non-Gaussian model [21, 22], where it was pro- 
posed that the fluctuations driving stock returns could be modeled by a 
'statistical feedback' process [23], namely 

dS = idSdt + aSdn (3) 

where 

dn = p{n)^duj. (4) 

In this equation, P corresponds to the probability distribution of fl, which si- 
multaneously evolves according to the corresponding nonlinear Fokker-Planck 
equation [24] 

dP _ 1 (9p2-<? 
'dt ~ 2 dfl^ ■ 
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The index q will be taken 3 > g > 1. In that case, Eq. (5) is known also 
as the fast diffusion equation [26]. Note also that in the above equation, the 
time t must be a-dimensional. In the following, the unit of time has been 
chosen to be one year. Correspondingly, a is also a-dimensional. 

Equation (5) can be solved exactly, leading, when the initial condition on 
P is a P{fl,t = 0) = d{fl), to a Student-t (or Tsallis [25]) distribution: 

p-^^i-^ + iQ-mm'it))-^' (6) 



with 
and 



l-q 



f3{t)^ctHi2-q){3-q)t)--^ (7) 



Z{t)^{{2-q){3-q)c,t)—^ (8) 
where the g-dependent constant Cg is given by 

2 



Cn 



/OO I 
{1 + {q - l)u'^)~^ du 
-OO 



(9) 



Eq. (6) recovers a Gaussian in the limit ? — > 1 while exhibiting power 
law tails for all g > 1. 

The statistical feedback term P can also be seen as a price-dependent 
volatility that captures the market sentiment. Intuitively, this means that 
if the market players observe unusually large deviations of fl (which - after 
removing a noise induced drift term which could equivalently have been ab- 
sorbed in the dynamics of Eq. (3) [21] - is essentially equal to the detrended 
and normalized log stock price) from its mean, then the effective volatility 
will be high because in such cases P{fi) is small, and the exponent q is larger 
than unity. Conversely, traders will react more moderately if Q is close to its 
more typical values. As a result, the model exhibits intermittent behaviour 
consistent with that observed in the effective volatility of markets - but see 
the discussion below. 

Option pricing based on the price dynamics elucidated above was solved 
in [21], and it was seen that those prices agreed very well with traded prices 
for instruments which have a symmetric underlying distribution, such as 
certain foreign exchange currency markets. However, the question of skew 
was not discussed in that paper and is instead the topic of the current work. 



5 



We extend the stock price model to include an effective volatility that is 
consistent with the leverage correlation effect, namely 

dS = fiSdt + a^o^-"5"d0 ( 10) 

with Q evolving according to Eq. (4). The parameter a introduces an asym- 
metric slcew into the distril:)ution of log stock returns. More precisely, when 
a < 1, the relative volatiUty can be seen to increase when S decreases, 
and vice- versa, an effect known as the leverage correlation (see [35]). For 
a — q — 1 the standard Black-Scholes model is recovered. For a = 1 but 
q > 1 the model reduces to that discussed in [21], while for q = 1 but general 
a it becomes the constant elasticity of variance (CEV) model of Cox and 
Ross [27]. 

There are two possible interpretations to model Eq. (10), and some limi- 
tations which we elucidate here. In the context of option pricing, the relevant 
question concerns the forward probability, estimated from now (t = 0), with 
the current price So corresponding to the reference price around which de- 
viations are measured. In this case, the fact that t = and S = Sq (or 
Q = 0) play a special role makes perfect sense. If, on the other hand, one 
wants to interpret Eq. (10) as a model for the real returns, then the choice 
of 5*0 = S{t = 0) as the reference price is somewhat arbitrary and therefore 
problematic. Still, this model produces returns which have many features 
consistent with real stock returns. This is illustrated in Figure 1 where we 
have plotted a time-series of simulated returns as well as, in Figure 2, a 
plot of the distribution of returns between times t and t + r, for different 
time lags r, averaged over all possible starting times t. Clearly, the model 
reproduces volatility clustering. Also, the returns distribution exhibits fat 
tails becoming Gaussian over larger time scales. (Note that this is not in 
contradiction with the fact that returns counted from t — have a Student-t 
distribution for all times t). The skew in the distribution is also apparent. 
For visual comparison, we show the same data for the SP500. The qualitative 
behaviour of the two data sets is very similar. However, in order to have a 
consistent model of real stock returns, one should allow the reference price 
to be itself time dependent. A possibility we are presently investigating is to 
write the statistical feedback term as P{Q — fi) (-"^"^^/^ , where Q is a moving 
average of past values of Q. (The current model corresponds to Cl — 0.) 
One can also easily alter the super- diffusive behaviour (as given by Eq(7)) 
of the distribution if a mean reversion of the fluctuations is included. Both 
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of these features make the model a more reahstic one of real returns (where 
distributions are non-Gaussian, yet volatilities are normally diffusive across 
all time-scales) . Incidentally, this super-diffusive volatihty scahng can alter- 
natively be modified by a redefinition of time (see also the discussion later 
in this paper). 

One point necessary to comment on regarding Eq. (10), is that depending 
on the value of a, stock prices may go negative. We do not allow this to 
happen: rather, we absorb the stock if the price hits zero. Basically, this 
means that the company has gone bankrupt, or that the stock is trading 
at such low prices (penny stocks) that it is practically worthless and we so 
not consider it part of the trading universe. It is interesting to look at the 
probability of bankruptcy V as a. function of a. This is shown in Figure 3. 
An analytical treatment [28], expected to be valid for rj — a{l — a) <^ 1, 
predicts that 

V^aT^T]^ (11) 

(where a is a computable constant), in good agreement with our data for 
small v. In fact, this point suggests a connection with credit risk markets, 
where the probability of default is an important quantity. One could imagine 
using a variant of our model to model default risk; conversely, one could 
imagine using data of probability of default to help choose the correct value 
of q; of a particular stock. The consequence of this default probabihty for 
option pricing will be discussed in Section 5. 

3 Fair Price of Options 

Our model contains only one source of randomness, cj, which is a standard 
Brownian noise. Therefore, the market is complete and usual hedging argu- 
ments are valid, as was shown in [21] for the case a — 1. It is however trivial 
to see that these results are also valid for arbitrary a. This means that we 
can immediately adopt many of the arguments usually associated with the 
Black-Scholcs log-normal world - even though wc here have a process with 
non-Gaussian statistics. In particular, it is possible to define a unique equiv- 
alent martingale measure Q to the process, related to the original measure 
F via the Radon-Nikodym derivative [29, 30], such that the discounted stock 
price G — Se~^^ is a martingale, where r corresponds to the risk-free rate. 
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12 3 4 

Time [Years] 

Figure 1: A typical path of the non-Gaussian model with q = 1.5 is shown. Here a = — 3 
and (7 = 30%. 

With respect to Q the process of Eq. (10) reads 

dS = rSdt + aSl~''S''dn (12) 

with 

dn = p'-^dz (13) 

where the noise z satisfies 

{dz{t')dz{t))Q = dtdt'S{t' - t) (14) 

Essentially then, we have replaced fi = r just as in the standard theory. Our 
task now is to solve for the option prices based on the model Eq. (12). 
The dynamics for G read 

dG = (7S'(5-"G'"e("~i)'^*cil] (15) 
= aS^-''G''p'-^e^''-''^''dz (16) 
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Figure 2: This plot exhibits the distribution of log returns {Y = InS) for the S&P500 
over different time lags r ranging from 1 day to 16, together with corresponding histograms 
obtained by our model with q = 1.5 and a = 1. 

We redefine time by introducing the new variable i such that 

2{a-l)rt _ 1 

Note that t —>■ t when {a — l)rt —>■ 0. With respect to t, the Brownian noise 
z satisfies 

{dzii')dzit))Q = e^^''-^^'-\dz{t')dz{t))Q (18) 

resulting in 

dG = aSl-'^G'^dn (19) 

d^{i) = p{Q{i))^dz{i) (20) 

Next, make the variable transformation 
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Figure 3: a) The probability of bankruptcy (default) P as a function of a for our model 
with q — 1.5, for times ranging from 0.25 to 2 years. The inset shows the probability of 
bankruptcy as a funtion of time for a stock with a = 0, compared to our prediction that 
V (X tor q = 1.5. These results are based on Monte-Carlo simulations with a = 30%, 
r = 6% and S'(O) = $50. b) The probability of bankruptcy is sensitive to the values of 
q and a as shown here, for 5(0) = $50, a = 0.2, r = 4% and T = 1. For T < 1 these 
probabilities are quickly depressed according to Eq(ll). 



(which simply becomes x = InG/^o for a = 1), so that 
dx{t) - 



a 



-a 



pl-q 

^ di + adn{i) 



2 1 + (1 - a)x{t) 
where the initial conditions read 

xo = x{0) = 

and where, at the real final time T, we have 



S{T) = e''^ (l + (1 - «)x(f )) 



(22) 



(23) 



(24) 
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Given a representation of the stock price with respect to the risk-neutral 
world in which the discounted price is a martingale, the fair value of a deriva- 
tive of the underlying stock can be calculated as the expectation of the payoff 
of the option. In this paper we shall focus on this approach. However we 
want to point out that one can equivalently obtain the option price by solv- 
ing the generalized Black- Scholes partial differential equation for the current 
problem, which can readily be written down following the lines of [21]. It 
reads 

where Pq evolves according to Eq. (6). In the limit q ^ 1 and a ^ 1, we 
recover the standard Black-Scholes differential equation. In the limit a — > 1 
we recover the case of [21]. 

In the following we proceed to study closed form option pricing formulas 
obtained via expectations. The most immediate and simplest path is to 
utilize our knowledge of the statistical properties of the random variable Q{t). 
By invoking approximations valid if cr^T <^ 1, which is certainly true for 
stock returns for reasonable maturities, we can obtain closed form solutions 
for European calls, much along the lines followed in [21]. However, we allude 
in Appendix B to a more complicated path which is based on mapping the 
process onto a higher dimensional process. 



4 Solutions via a Generalized Feynman-Kac 
Approach 

If a^T <^ 1, which is valid even for relatively high volatility stocks, (e.g. 
typical volatilities of 10% to 30% yield a^T values of .01 to .09 for T = 1 
year), then we can insert the approximation x{t) ~ a^}{t) into the right hand 
side of Eq. (22) yielding, 

x{f) = an{f) - r (26) 

plus order corrections. For a = 1 and general q, the problem reduces to 
that discussed in [21], which we shall revisit with a slightly different evalua- 
tion technique below. 
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The integral in Eq. (26) contains terms of the type 



rF{n{t))dt, (27) 
Jo 



in other words integrals of a function F of the random path 0,{t), conditioned 
on ending at a particular value f2(T) = Q^- One way to evaluate the integral, 
which is the approach taken in [21], is to invoke the following (which is 
only valid approximately [31]): Replace the actual paths Q{t) with other 
paths, such that the ensemble average of the two sets of paths are the same. 
Such paths can be obtained by exploiting the scaling properties of the time- 
dependent probability distribution P{Q{t)). We know that the variable Q{t) 
at time t is distributed according to Eq. (6) above. Due to scaling, we know 
that we can map Q{t) onto the terminal value Qt in the following way: 



fl{t) 



where ^{T) is distributed according to Eq. (6) evaluated at time T. There- 
fore, as we illustrate in Appendix A, the ensemble statistics of this replace- 
ment path and the original path are equivalent. This approach was shown 
in [21] to be quite precise numerically. However, as we discuss next, there is 
another more exact approximation which can be used to evaluate Eq. (27). 
This entails solving a Feynman-Kac equation exactly up to first order in cr^. 
Define the quantity 

Wi^T, T)^y^Wi exp{e C F{n{t))dt} (29) 

which is the expectation calculated as the sum over all paths i ending at Jlr 
at time T, of the exponential of quantity which we wish to calculate. The 
coefficients Wi denote the weights of the paths, and e is a small parameter. 
The function F can be a general one of ^l{t), but in our current problem we 
shall be only interested in polynomial expansion up to the second order in 
Q. We can write 

WiQr, T)^Y.^i exp{e / (E fjm{ty)dt} (30) 
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In the case where the paths evolve according to the nonhnear Fokker-Planck 
equation, one can establish a generalized Feynman-Kac equation for W{flT, T) 
that reads: 



W ~ 2 



eE/.(^)^T (31) 
j 

where Wo — Pq{flT,T) represents the solution to the fast diffusion problem, 
corresponding to e = 0. It is relatively straightforward to insert the Ansatz 

W^Wo[l + e{go{T) + g,{T)nT + 92{T)nl)] (32) 

into Eq. (31) and obtain an exact set of differential equations for the time 
dependent coefficients gj, given in Appendix A (Eq. (52) - Eq. (54)). 

First, let us illustrate this approach on the example of a = 1 which cor- 
responds to the problem discussed in [21]. The expression for S{T) becomes 
(from Eq. (24) and Eq. (26) in the hmit a ^ 1) 

Sm = S(Q) expirT + a^r - ^ I Zitf-ni + iq- l)(3(m(t?)dt} (33) 

2 Jo 

Using the Feynman-Kac formula to evaluate 

rT 



J 

Jo 



Z{tf-''(3{t)VL{tfdt (34) 



(see Appendix A, Eq.(56) with e = (T^/2 and /i2 = ^/3), we obtain an 
expression of the form: 



( r Zitf-^mm^dt) = go{T) + g^{T)nl (35) 
Jo 
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{gi is zero by symmetry in this case). For the stock price, we thus obtain 

^2 r 3 _ 5 



S{T) = S{0) exp{rT + a^T - ^ 



l{T)^ + {q-l){go{T)+g,{T)^l'^) 



} 

(36) 

with 7 and g^ are given by Eq. (57) - Eq. (59). The above expres- 
sion is exact to order o"^; to that order, the path to path fluctuations of 

/o^ Z{tf-'il3{t)Vt{tfdt conditioned on a given value of VLt can be neglected. 

The coefficients gj found here are slightly different than those in [21], 
where instead of Eq. (34), the approximation Eq. (28) was used. However, 
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the results are numerically very close, so that in practice either evaluation 
method can be used; both give option price which match well to Monte-Carlo 
simulations. In fact, why this is so becomes clear from the plots shown in 
Appendix A, where Monte-Carlo simulations of the quantity Eq. (34) are 
compared with the Feynman-Kac approximation and the one of [21]. In 
the general case where a is not restricted to unity, the challenge now is to 
evaluate Eq. (26) which can be rewritten as: 



We shall focus our discussion on the last term. A solution can be found 
using a two-step approach: First of all, we make the assumption that the 
expression, when integrated to abitrary time m, can be well-described by a 
Padc approximation, namely a ratio of polynomials, which has the correct 
asymptotic behaviour for Q — > and Q ^ oo, resulting in 



Secondly, the coefficients can be determined by equating the Pade expansion 
in the small and large Qj, limits with the corresponding Feynman-Kac 
expectations of simple polynomials. The details of this calculation are in 
Appendix A, and yields values of the coefficients A, 5, C and D as given 
in Eq. (68) below. Again note that even in this general case, the naive 
approximation of Eq. (28) yields a slightly different result, yet numerically 
the two are practically indistinguishable. This can again be understood by 
the plots shown in Appendix A, comparing the two approximations with 
Monte-Carlo simulations. Note also that the Pade approximation of order 
2:1 which we use in Eq. (37) is already extremely close to the Monte-Carlo 
results (see Appendix A, Figure 12b); yet if desired, convergence can easily 
be further improved simply by including higher order terms. 

5 European Call 

The final expression for the stock price at time T with respect to the mar- 
tingale noise is thus given by Eq. (24) together with Eq. (38). With this 




(37) 



x{f) = a^f 



aa^A{T) + B{T)nf + C{T)n^^ 

2 i + D{f)nf 



(38) 
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result, it is straightforward to obtain a general expression for the price of a 
European call in the current framework. The price of a European claim c 
can formally be written as 

c = e-^^{h{ST))Q (39) 

where h is the payoff of the option. For a European call the payoff is h — 
max(5'T — K,0), allowing us to evaluate the fair call price as 

c^e~^^{ST-K))^ (40) 

where D represents the domain of non-zero payoff, namely 

St> K (41) 

In our current framework there must be an additional constraint on the paths 
expiring in the money, namely that they never crossed S{t) = at any t < T. 
The probability p = P{S{T) > K \ S{t) < 0) that any path would cross and 
then return to expire greater than K is intuitively the order of the probability 
of default squared. More exactly though, it can be computed in the limit 
where r] = a{l — a) is small, corresponding to small default probabilities V. 
Using most probable path methods [28], we find 

P«^xexp(-^), (42) 

where 6 is a positive constant of order unity. For practical purposes, this 
difference is extremely small when compared to the accuracy of our solution, 
and can be neglected in most cases. 

The condition St — K implies a quadratic equation for which can 
easily be solved, resulting in two roots di and ^2 given in Appendix A, Eq. 
(69). These define a domain of integration for which the inequality St > K 
is satisfied. It can be seen that di always corresponds to the relevant lower 
root, while the upper bound is dmax = '^2 if d2 > di, and dmax = 00 otherwise. 
(Note that the strictness of this upper bound is however in principle irrelevant 
because of the smallness of the probability distribution in that region, and 
also because our approximation scheme breaks down in that region). 

We are now ready to state one of our main results, namely the closed- 
form expression for the price of a European call option within this framework. 
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Following Eq. (40), it is given as 



c ^ So 



/ {l + {l-a)x{T)) — Pg{Qf)dnf 



f 

Jd 



Pg{Qf)dflf 



(43) 



with Xf a function of flf given by Eq. (38) evaluated at m = T, T is given 
by Eq; (17) with t — T, and Pq(Q^) is the probability of the final value of il, 
given by: 



Note that it is at this point that we are in fact overcounting those paths 
which crossed zero yet expired above K. To account for this, Pq should be 
replaced with Pq —p of Eq. (42). However, as mentioned above, p can readily 
be dropped as it leads to negligibly small contributions, and the above result 
for the option price can be seen as exact to order a^. 



As already mentioned, several special cases are recovered for certain values 
of q and a. In practice, occasions could arise in which it might be useful to 
work with these simpler solutions. Therefore, we spend a few lines discussing 
them here. 



The case oi q — 1 and general a corresponds to the CEV model of Cox and 
Ross [27] . One obtains for the call price 




(44) 



6 Special Cases 



6.1 CEV model q=l 



/ {l + {l-a)x{T)) — p{zf)dzf-e-''^K 



C— Si 







Jdl 



/ P{zf)dzf (45) 



with notation zt = Qq=i{t) and where 



, ^ , a; n - 
x{T) - --a^T + r 



(46) 
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with 



and 



r = 1 + ^Min^^f (47) 



Of course, Zt is a Gaussian variable so 

These equations look different to the solutions presented in [27], which are 
expressed in terms of Bessel functions. Numerically, however, our solution 
and theirs are the same. The reason that they look different is - apart from 
the small approximation we have made - is that we have made a point 
of explicitly averaging with respect to the distribution of the random noise 
variable z rather than with respect to the distribution of the stock price itself. 
The reason for this is because as we saw above in the general case, we know 
something about the distribution of the noise Q, whereas the distribution of 
S itself is more complicated - but see Appendix B. 



6.2 Non- Gaussian Additive or q-Normal Model a = 

The CEV model with a = is known as the normal model. Here wc present 
the solution of the corresponding model for general q, which wc term the 
q-normal model. In this particular case, it is easy to solve for the option 
price in terms of the distribution of S, since the noise is additive. The call 
price reads: 



-rT 



Z{T)a Jk 



S 



-imts- 



l + (9-l)^(S-5oe'-T 



a 



q-l 



dS 



Z{T)a Jk 



dS 



(50) 
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6.3 Non-Gaussian Multiplicative a = 1 



Finally, the case studied in [21] is recovered for general q, and o = l. For 
completeness we cite this result, albeit using the more precise Feynman-Kac 
evaluation of the path integral which resulted in Eq. (36): 

c = 5o/^^exp|al]^-^[7(T)^-(l-g)(<7o(T) + <72(T)(^|)]| 

Pg{nT)dnT - e-^'^K / Pg{nT)dnT (si) 

with Pq as in Eq. (44) and di is the smallest root of the equation Eq. (36) 
such that S{flT) = K. 



7 Numerical Results 

Now that we have pricing formulas for a non-Gaussian model with skew, it 
is interesting to look at the option prices which result from the model, and 
to see how they compare with a standard Black-Scholes formalism. Using 
Eq. (43) with q = 1.5, a = -1.5, a = 30%, r = 6% and = $50, we 
calculated call option prices as a function of different strikes K and times 
to expiration T. We then backed out imphed volatihties, i.e. the value of 
a which would have been needed in conjunction with the standard Black- 
Scholes model {q = a = 1) to reproduce the same option prices. A plot 
of those implied volatilities as a function of K and T constitutes the skew 
surface which we show in Figure 4. The general shape of this surface is very 
similar to what is observed on the marketplace. For small T, the profile 
across strikes is smile-like, becoming more and more of a downward sloping 
smirk as T increases. In Figure 5 we show a similar plot for q = 1 (i.e. the 
CEV model). Note that this skew surface does not capture the shape one 
observes empirically. This goes to show that both fat tails {q > 1) and skew 
(a < 1) are necessary for a realistic description. 

It is also interesting to look at the variation of the option price with 
respect to the new parameter a (the variations to q, denoted by Upsilon T, 
were studied in [21]). We denote this new 'Greek' by the Hebrew letter Aleph 
K, such that K = dc/da. Figure 6 shows K as a function of stock price 5*0, 
using K = $50, T — 0.5 years, a — 30%, and r = 6%, for q — 1.5 and q — 1. 
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Figure 4: A plot of the skew surface, i.e. Black-Scholes implied volatilities across strikes 
(K) and time to expiration (T), backed out from our non-Gaussian model with q — 1.5 
and a ~ —1.5. Other parameters were Sq — $50, r — 6% and a — 30%. 

The asymmetric nature of the sensitivity to variations in Aleph is apparent 
for q = 1.5. Another way of depicting much the same information is shown 
in Figure 7, where volatihties imphed by comparing standard Black-Scholes 
with our model Eq. (43) for different values of a are shown. In all cases 
we used q = 1.5, Sq = 50, T = 0.5 years, cr = 30% and r = 6%. One can 
clearly see that the implied volatility curve is like a smile for a = 1 (no skew) 
becoming more and more asymmetric about = 5*0 as a decreases. 

In all of these results, we used the closed form pricing formula, Eq. (43), 
to generate option prices. However, since we used some approximations along 
the way, it is a good check to see how the closed form price compares with that 
obtained from pricing via Monte-Carlo simulations of the process. Indeed, 
we saw that the two values are indistinguishable within the limits of accuracy 
of the Monte-Carlo simulations. 



19 



Figure 5: The skew surface implied from a CEV model {q = 1) with a — —1.5. Other 
parameters were 5*0 — $50, r — 6% and a = 30%. The kink at short times may just be a 
numerical artefact. 

8 Empirical Results 

To illustrate the qualitative agreement between our model and market data, 
we show in Figure 8 a plot of the empirical skew surface for OEX options, as 
well as the skew surface backed out of our model with q = 1.5 and a = —1.2. 
We have not at all tried to calibrate the model to match the empirical data 
in any way, the plot is only intending to show that our model produces a 
surface with similar features to the empirical one across several time scales, 
with just one set of parameters a, a and q. In other words while with respect 
to the Black-Scholes model, the entire skew surface of Figure 8b is needed 
to describe the option prices, with respect to the q = 1.5 and a = —1.2 
model the skew surface reduces to a single point. Consequently, we have 
the hope that real market smiles and skews might be captured by our model 
with parameters q, a and a varying only slightly across strikes and times to 
expiration. 

Clearly, to test this with any statistical significance requires a large study 
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Figure 6: Aleph the Greek, H: The partial derivative of the call price with respect to a 
is plotted as a function of the stock price 5(0), for q — 1.5 and q = 1.. The asymmetric 
nature is apparent for q = 1.5. We used K = $50, r = 6%, a — 30% and T — 0.5 years. 

on many options, which we leave for a future work. However, we do present, 
again for the purpose of illustration, results based on an analysis of one set 
of options, namely call options on Microsoft (MSFT) traded on November 
19, 2003. A popular methodology that market makers and traders follow 
is to vary the parameters of whatever model they are using such that the 
theoretical smile matches the market at each time to expiration T. We could 
also follow such an approach: vary q, a and a for each value of T such that 
the smiles and skews are reproduced. But if the model is good in the sense 
that the parameters do not change much over time, then perhaps the most 
parsimonious treatment would be to choose one set of those parameters such 
that the entire skew surface across both strikes and expiration times is well- 
fit. While interesting, both of these approaches would be ways of implying 
the model parameters from the options data. 

An alternative methodology which is well-suited for our current approach, 
is to instead try to relate at least one or two of the model parameters to prop- 
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Figure 7: Implied volatilities from the q = 1.5 model as a function of a, with Sq = 
$50, r = 6%, a = 30%, and T = 0.5 years. As a decreases, the skew increases. Intuitively 
this make sense: a contributes to larger negative tails in the distribution of the underlying, 
and therefore there will be higher probability of expiring out of the money relative to a 
Black-Scholes process with lognormal noise. For extreme out of the money values however, 
the noise in the fat tails can again increase the probability that an out-of-the money option 
can expire in the money, so we expect the smile to increase again in this regime. 

erties of the underlying asset, and then use these to predict market smiles. 
In particular, the parameter q can be readily determined from the empirical 
distribution of the returns of the underlying asset. It has been found in pre- 
vious studies [33, 21, 34] that a value of g ~ 1.4 captures well the distribution 
of daily stock returns, so this value could be adopted in the option pricing 
formula. The parameter a could in principle also be determined from the 
distribution of underlying returns, or from measuring a leverage correlation 
function [35, 10]. Alternatively, it could perhaps be determined (as men- 
tioned earlier in this paper) form empirical probabilities of default. For the 
present study we fix only q from the underlying distribution, and imply a 
and a from the market smiles. We shall then look at how good the implied 
smiles fit the data, in conjunction with how the implied parameters vary with 
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Figure 8: A purely qualitative comparison between the empirical skew surface of a set 
of OEX options on S&PlOO futures traded on June 6, 2001, and the implied skew surface 
from our model with q = 1.5, a = —1.2, a = 30%, r = 4.5% and Sq = 660. We have not 
tried to calibrate to the OEX data in any way, we simply wish to show that the general 
behaviour of the surfaces across strikes and time to expiration is similar. From top to 
bottom: T = 0.03, T = 0.12, T = 0.20, T = 0.29 and T = 0.55. 



T. If they are somewhat stable then we can conclude that the model is quite 
good. 

In Figure 9 we show the results for our model. As T ranges from the 
order of a month to a year, we choose a and a such that the best fit in 
terms of least square pricing error is obtained between model and empirical 
call prices, for fixed q. The plots show actual implied volatility and the 
imphed volatility of our model. It is clear to see that the q — lA model 
provides a good fit of observed smiles at each time to expiration. However 
for the longest maturity T — 1.17 one sees that the away-from-the money 
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Figure 9: A quantitative comparison between the empirical skews obtained from a set of 
MSFT options traded on November 19, 2003, and our model with q= 1.4, which well-fits 
the returns distribution of the underlying stock. We varied a and the volatility parameter 
a of the q = 1.4 model for each time to expiration T (see Table 1). We used r = 4.5% and 
Sq = $25.55. AIV stands for the actual average (of put and call) implied volatility. 

strikes are slightly over-valued by our model. The valid question remains as 
to whether this discrepancy is a true market mispricing or an artifact of the 
model, which does not lead to log-normal statistics for large T. (Note that 
we can easily match the market exactly even at T = 1.17 simply by reducing 
q appropriately, but this is not what we are attempting to do here.) 

In Table 1 we show how the parameters of the model vary as a function 
of T. The parameter a goes from 0.1 at early times to fluctuate around 0.2. 
Note that this order of magnitude of a, when related back to the distribution 
of stock returns, is entirely consistent with empirical observations of the 
leverage effect [35]. The volatility parameter a exhibits a negative term 
structure, ranging from 32% to 25%. Such a negative term structure has 
been consistently observed in our empirical studies (for example it is seen in 
ref. [21] where options on FX futures are analyzed). The reason for such 
an effect is, as mentioned earlier, that the g-model considered here predicts 
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T [Years] 


0.082 


0.159 


0.41 


1.17 


a[%] 


32 


31 


27 


25 


a 


0.1 


0.2 


0.3 


0.2 



Table 1: Variation of the fitted parameters a and a with time to expiration T, for q= 1.4. 

an anomalous growth of the volatility with time, as T^/^'^~i\ instead of the 
empirically observed standard diffusive scaling, a/T. A way to correct for 
this is to allow the volatility a to be maturity dependent, which amounts 
to a mere redefinition of time in the nonlinear Fokker-Planck equation that 
defines the model. For g = 1.4, a should scale as T~^/^ to reproduce the 
correct time dependence of the width of the return distribution with time. 
This scaling appears indeed to be satisfied as shown in Figure 10, where 
we plot a (Table 1) of the MSFT data against time on a log-log scale. To 
illustrate the regularity of this temporal behaviour, we also show in Figure 
10 the corresponding plot for the a parameter of the options on FX futures 
initially shown in [21]. Based on these results, we see that the total variation 
of the model parameters a and the maturity rescaled a are only very slight. 
Furthermore, by construction q is kept constant. 

Another comment, which should be taken loosely since we have not done 
a systematic study of this point, concerns the probabilities of default implied 
by the parameters of this MSFT example, when reinserted in our stock price 
model. Assuming that the real drift of MSFT is the risk free rate (which is 
certainly an underestimate), one obtains < 0.005% probability of bankruptcy 
for T = .082, 0.01% for T = 0.159, 0.07% for T = 0.41 and 0.35% for T = 
1.17. These estimates were obtained from Monte Carlo simulations, and do 
not seem unreasonable. However, the above probability of default at 1 year is 
probably higher that typical credit risk ratings of MSFT. This overestimation 
of the probability of default at longer timescales is intimately related to the 
smile shown in Figure 9, for T — 1.17: the away-from-the money strikes are 
slightly over-valued by our model as T increases. As mentioned earher, this 
is because we choose to keep q = 1.4 fixed at all timescales. If we instead 
decided to vary q to fit the smile exactly, it would be closer to g = 1 and 
the default probabilities would drop much closer to again (for example, see 
Figure 3 b). Also note that in reality, default probabilities might depend on 
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Figure 10: Deterministic term-structure in the volatility parameter a of the g = 1.4 
model, shown (top) for the MSFT options example of Figure 9 and Table 1; shown (bot- 
tom) for the JY Futures options example discussed in [21]. 

the real rate of return of the stock, which for MSFT is substantially higher 
than the risk free rate. 

It could be interesting to place our results in the context of a comparison 
with a popular stochastic volatihty model (namely the SABR model, which 
for the sake of this discussion is briefly summarized in Appendix C). The 
fit of the SABR model to the MSFT options discussed above was performed 
by an external source [32] so we only briefly report the results here. Three 
parameters were varied at each T, and consequently the smiles could be fit 
much as in our example. The volatility parameter a varied with a slightly 
positive term structure around a value of 29%. However the parameters p 
(related to the correlations between stock fiuctuations and volatility correla- 
tions) and A (which describes the volatility of the volatility) did vary quite 
a bit. We saw that p increased from 0.1 to 0.6. Perhaps more significantly, 
A went from around 9.15 to 1.14 as maturity increased from 0.082 to 1.17. 
This decrease makes sense because at large times, the stochastic volatihty in 
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the SABR model is unbounded; consequently the implied option smile would 
be way too pronounced if A did not suppress this feature, forcing the model 
into the log-normal limit as T increases. On the whole, we saw that while the 
SABR model certainly fit the market smiles well, the number of parameters 
(3 if (3 is fixed) and their variation is larger than in the non-Gaussian q — lA 
model (with 2 free parameters). 

Consequently, we can conclude (at least for the options studied in this 
example), that our non-Gaussian model with fixed q indeed yields a rela- 
tively parsimonious description of empirically observed option prices, cap- 
turing both the smile and the skew with a few seemingly robust parameters. 

9 Conclusions 

In this paper we have extended the non-Gaussian option pricing theory of 
[21] (which recovers the standard Black-Scholes case in the hmit q = 1) 
to include asymmetries in the underlying process. These were introduced 
so as to incorporate the leverage correlation effect [10] in a fashion such 
that the CEV model of Cox and Ross [27] is recovered in appropriate limits 
{q = a = 1). A theoretical treatment of the problem is possible much along 
standard lines of mathematical finance. Using the fact that the volatility 
of the process is a deterministic function of the stock value, the zero-risk 
property of the Black-Scholes hold, and one can set up a generalized Black- 
Scholes PDE as well as define a unique martingale measure allowing us to 
evaluate option prices via risk-neutral expectations. We have introduced 
a generalized Feynman-Kac formula for the family of stochastic processes 
involved in our model, namely statistical feedback equations of the type [23] 
which evolve according to a nonlinear Fokker-Planck equation [24], [26]. This 
formula in conjunction with Fade expansions could then be used to evaluate 
closed-form option pricing equations for European call options. 

Numerical results allow us to back out Black-Scholes implied volatilities. 
A plot of those across strikes K and time to expiration T constitute a skew 
surface. We found that the skew surfaces resulting from our model exhibit 
many properties seen in real markets. In particular the shape of the surface 
tends to go from a pronounced smile (across strikes) to a sloping line as T 
increases. A comparison of the model to a restricted set of MSFT options 
was discussed. We found that the entire skew surface could be well-explained 
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with a fixed value of g = 1.4, whicfi also fits to the distribution of the under- 
lying. The remaining parameters a and o vary only slightly, especially after 
we factored out a deterministic term structure in the volatility parameter. 
While the philosophy of many market participants is in general to tune the 
parameters of their models to fit market smiles, we try to relate our param- 
eter to the underlying distribution in the hope that we can attain a more 
parsimonious (and therefore more stable) description of both underlying and 
option. (Along this hue of thought, see also [36] and refs. therein). Indeed, 
our results - though not yet statistically significant - do indicate that it might 
be possible to explain the entire skew surface with one set of constant (or 
slowly varying) parameters. We hope to strengthen this statement through 
more empirical studies. 

Acknowledgements: We wish to thank Jeremy Evnine, Roberto Os- 
orio and Peter Carr for interesting and useful comments. Benoit Pochart is 
acknowledged for careful reading of the manuscript. 



10 Appendix A: Feynman-Kac Expectations 
and Fade Coefficients 

The coefficients of the Feynman-Kac Ansatz Eq. (32) when inserted into Eq. 
(31) must satisfy 

^ = Z{ut-^g^{u) (52) 

^ = 2{q-2)Z{ut-^i5{u)9x^hM (53) 
au 

^ = {hq-^)Z{uf-'^ii{u)g^{u)-rh^{u) (54) 
au 

From these equations follows that for path integrals of type exp{ /q" h\ {t)fl{t)dt} 
(with h2 — 0), then only the gi coefficient is relevant and the expected value 
of the integral yields the approximation 



r hAt)Vtit)dt ^ gAu)Vtu (55) 

^0 
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Similarly, if we are only looking at path integrals of the type expl/g" h2{t)VL{tY di} 
(with hi = 0), coefficients go and g2 are coupled together, so that the 
Feynman-Kac approximation of the integral implies 



h2m{tf ^ go{u) + g2{u)nl (56) 



In the option pricing problem without skew, namely the case a = 1 of 
Eq. (33), we have hi — and /i2 = Integration yields 



92[U) 



2(9 - 5q) 



9-5g 



(57) 
(58) 



with 



^{u) = ((3 - g)(2 - q)c,)—^u—^ (59) 

For the general case of option pricing including skew (general a as in Eq. 
(37)), we expand the Fade Ansatz of Eq. (38) for both small and large Q, 
and equate with the same expansions of the actual quantity we are interested 
in, which we then evaluate using Feynman-Kac expectations. In the small fl 
case: 



A + {B-AD)n^+{C-BD + AD'^)fll 



(60) 



= / z(ty-''\i-rin(t) + (i-q)(3{t)n{tf 

Jo L 

= r \z{ty-i - 7]hi{t)n{t) + (1 - q)h2{t)n{tf 

Jo L 



dt 



- + (1 - q){go{u) + g2{u)nl) 



with 

77 = a(l - a) (61) 

and 7 as in Eq. (59). 

The coefficients g-j arc calculated from Eq. (52)-Eq. (54) with hi = Z^~^ 
and /i2 = PZ^^'^ and result in: 
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3-g 

2(9 - 5q) 



(62) 
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Si = 7(»)^^ (63) 



Similarly, the large Q expansion yields 



= - — - I h{t)n{t)dt (65) 



D r] Jo 

= - — -gi{u)flu (66) 
V 

with gi calculated from Eq. (53) using hi — l3{t)Z{ty~'^, yielding 

Through standard coefficient comparison, we have enough information to 
solve for the Pade coefficients A,B,C and D, resulting in 

^ " '-7 (68) 



A = 


^o(g-i) + 


B = 


AD - r]gi 


C = 


(q-lf-D 
V 


D = 


92{q - 1) 





In the case of a general skew, the condition St — K yields a quadratic 
equation with the roots 



with 



2M ^ ' 



1 — a 2 

2 

M = C^-aD (71) 
R = <^''"/"°'"°-V^^ (72) 
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where A, B,C and D are evaluated from Eq. (68) at the time u = T with f 
ofEq. (17). 

These results are all based on evaluating the terms of type Eq. (27) using 
the Feynman-Kac Ansatz Eq. (32). We want to compare in details these 
results to the ones obtained in [21] using a different evaluation technique, 
namely replacing the path Q[t) by another path defined such that the en- 
semble distribution at each time would be the same as for that of the original 
paths, namely Q{t) = ^ [3{T) / [3{i)VLT as in Eq. (28). The numerical results 
obtained by that approximations are extremely close to those obtained by the 
current evaluation method, and by Monte-Carlo simulations. To elucidate 
this point we show in Figure 11a a plot of the quantity 

r l3it)Z^-mitfdt (73) 
Jo 

versus Q,t evaluated for q — 1.5 and T — 0.5 via i) Monte-Carlo simulations, 
ii) as in Eq. (28) and iii) using the Feynman-Kac evaluation Eq. (35) with 
Eq. (57) and Eq. (58). It is clear that the approach Eq. (28) slightly under- 
values the expectation of the true paths for small i^(t), and over- values for 
large Vt{t), in such a way that the average value over all Vt{t) yields the correct 
result. The Feynman-Kac approximation, on the other hand, is expected to 
be exact in this case. 

In the more general case with a skew, we need to evaluate an expression 
of form 

Again this can be done for q — 1.5, a — 0.5 and T — 0.5 using the i) Monte- 
Carlo simulations, ii) Eq. (28), iii) Feynman-Kac equation together with a 
Fade expansion. The results are shown in Figure lib. Similar behaviour as 
that of Figure 11a is exhibited: the effective path approximation of Eq. (28) 
under- values and over- values the true result in such a way that when averaged 
over all a good approximation is obtained, whereas the Feynman-Kac 
approach is uniformly better. 

The close agreement between both approximation techniques allows us in 
practice to use either one for option price evaluation. The numerical results 
are extremely close. Nevertheless, the Feynman-Kac approach should be 
preferred as the more exact one. 
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Figure 11: Evaluation of path dependent integrals as a function of the terminal value r^T, 
using i) Monte-Carlo simulations, ii) the approximation of Eq. (28), and iii) Feynman-Kac 
techniques, a) The quantity Eq. (73) relevant for q = 1.5 and a = 1.0 (no skew) and b) 
the quantity Eq. (74) relevant for the general skew case, here with q ~ 1.5 and a — 0.5. 
In both cases cr = 30%. 

11 Appendix B: Exact Solutions via Hyper- 
Geometric Functions - A Proposal 

As an addition to the theoretical part of our paper we would like to briefly 
discuss a possible alternative path to solving for the option prices of the 
non-Gaussian model with skew Eq. (10). The solutions will involve the fact 
that one can map the transformed problem of Eq. (22) onto a free-particle 
problem in higher dimension as described below. 

In the standard case q = 1 the CEV model of Cox and Ross admits an 
explicit solution in terms of Bessel functions for arbitrary values of a. In 
this appendix, we show how this result can be obtained by mapping to CEV 
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process to a standard Brownian motion in higher dimensions, and how this 
method generahzes, although incompletely, to the case q > 1. 

Starting from Eq (22), the change of variable y — 1 + {1 — a)x leads to: 

dy = rjdn - — -rj'P'-^y-'dt (75) 

z(l — a) 

with 77 = (t(1 — a), in this appendix we drop the hat on the time variable. 
To order 77^, this corresponds to a Fokker-Planck equation of form 

^ - a—^^ + (76) 
dt dy y dy"^ 

where the time has been rescaled by rf and a = a/{l — a). 
Now, inserting the Ansatz P — f{y,t)^o{y) one obtains 

+ - i)r-\^)'' + I'tr'r^ + (77) 



with u = 2 — q. 

Grouping together all terms which do not contain derivatives of /, and 
imposing that the coefficient vanishes leads to an ordinary differential equa- 
tion for $0^ 

which is solved by 

$0 = (79) 

provided A(i/, a) satisfies the following quadratic equation: 

ai/^A^ + uX{a - 1) - a = 0. (80) 

The correct root is the one that vanishes when a — > 0. All remaining terms 
can be regrouped as 
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where 

A - ^ + lul-^ (82) 
dy^ y dy 

is the radial Laplacian operator in d (fictitious) dimensions, where d is given 
by: d—1 = a+2Xi'. The y-'^^^^i) term can finally be eliminated by introducing 
a new coordinate g — g{y) defined as: 

l_A(,.-l)/2 

«= l-A(.-l)/2 - 

In terms of this new coordinate, / obeys a non linear radial diff'usion equation 
in dimensions d' 

% = A,,ng,t) (84) 
where the efi^ective dimension d' is finally given by: 

2-A(.-l) ^''^ 



Let us first focus on the case q = v = 1, where the diffusion equation is 
linear. The initial condition on x, xq = 0, translates into an initial condition 
on f{g) which is a (5-function over the hyper-sphere in dimension d' — d, of 
radius go — g{0)- The solution of the radial diffusion equation in d dimension 
for an isotropic initial condition is obviously constructed as the superposition 
of point source solutions of the standard d dimensional diffusion equation (i.e. 
a d dimensional Gaussian), averaged over the position of the starting points, 
here sitting on the hyper-sphere »So of radius go- More explicitly, introducing 
d dimensional vectors g, one has: 

= (86) 
Introducing the angle 9 between g and one finds: 

^^'^ ^ (4^^ [-^T ) Jo ^ ] ■ ^''^ 
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Using the following identity: 

Jo ^ = ^^^^^ fcj ''^z-^^^^' 

(88) 

where / is the Bessel function, we finally recover the solution of Cox and 
Ross (see also [37]). 

The case q > 1, u < 1 leads to the so-called 'fast' diffusion equation in 
d' dimensions. In this case, an explicit point source solution can be easily 
constructed for an arbitrary position of the point source, and is similar to the 
d' = 1, Student-like solution discussed in the main text. Unfortunately, the 
general solution for an arbitrary distribution of point sources can no longer 
be constructed since the equation is non-linear. The case where these points 
are on an hyper-sphere is, to the best of our knowledge, unknown, although 
approximate solutions could perhaps be constructed both for short and long 
times. We leave the investigation of this path for future work. 

12 Appendix C: The SABR Model 

The SABR model [6] is a stochastic volatility model of the following form 

dS = S^aduui (89) 

da = Xadu;2 (90) 

with 

< diiJiduj'i >— pdt (91) 

and 

a(0) = (7 (92) 

It contains four parameters, (5 and p both contribute to the skew {(5 is in fact 
what we call a in the present paper), while A is related to the curvature of 
the smile and a is the volatility parameter. Because the log-volatility follows 
a purely diffusive process it can become arbitrarily large at large times. In 
our study we assume /3 held fixed, seeing as varying p can already change 
the slope of the skew curve. 
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